#include <stdio.h>
#include <pthread.h>
#include <stdlib.h>
#include <math.h>

#include "gauss_elim.h"

struct drand48_data g_buffer;
double **gMatrixA;
double gVectorB[SIZE];
double gVectorX[SIZE]={0, };

int getMaxLoc(int nStartRow, int nCol);
void showMatrix();

/**
  * Initialize program 
  * - Generate matrix A & array B
  */
int init()
{
	int i, j;
	/* Memory allocation */
	gMatrixA = (double**)malloc(sizeof(double *)*SIZE);
	for(i=0; i<SIZE; ++i)
		gMatrixA[i] = (double*)malloc(sizeof(double)*SIZE);

	/* TODO: Need to parallize */
	srand48_r(time(NULL), &g_buffer);
	for(i=0; i<SIZE; ++i){
		for(j=0; j<SIZE; ++j){
			drand48_r(&g_buffer, &gMatrixA[i][j]);
		}
		drand48_r(&g_buffer, &gVectorB[i]);
	}
	
	return 0;
}


/** 
  * Finilize this program
  */
int fini()
{
	int i;
	for(i=0; i<SIZE; ++i)
		free(gMatrixA[i]);
	free(gMatrixA);

	return 0;
}

int gaussian_elim()
{
	int i, j, k, nMaxLoc;
	double *pTmp, dCoeffi;
	for(i=0; i<SIZE; ++i){
	/* TODO check whether this part really necessary or not */
#if 0
		nMaxLoc = getMaxLoc(i, i);
		pTmp = gMatrixA[i];
		gMatrixA[i] = gMatrixA[nMaxLoc];
		gMatrixA[nMaxLoc] = pTmp;
#endif	
		/* TODO Need to parallize */
		for(j=i+1; j<SIZE; ++j){
			dCoeffi = gMatrixA[j][i]/gMatrixA[i][i];
			for(k=0; k<SIZE; ++k) 
				gMatrixA[j][k] = gMatrixA[j][k]-dCoeffi*gMatrixA[i][k];
			gVectorB[j] -= dCoeffi*gVectorB[i];
		}
		// showMatrix();
	}
	return 0;
}

int getMaxLoc(int nStartRow, int nCol)
{
	double dMax=-1.0;
	int nMaxRow=-1;
	int i;
	for(i=nStartRow; i<SIZE; ++i){
		if(dMax<gMatrixA[i][nCol]){
			dMax = gMatrixA[i][nCol];
			nMaxRow = i;
		}
	}
	return nMaxRow;
}

int backSubstitution()
{
	int i, j;
	double dTmp;
	for(i=SIZE-1; i>=0; --i){
		dTmp = gVectorB[i];
		for(j=i+1; j<SIZE; ++j){
			dTmp -= gVectorX[j]*gMatrixA[i][j];
		}
		gVectorX[i] = dTmp/gMatrixA[i][i];
	}
	return 0;
}

int getL2Norm(double* pL2Norm)
{
	double dL2Norm=0;
	double dTmp;
	int i, j;
	for(i=0; i<SIZE; ++i){
		dTmp = 0;
		for(j=0; j<SIZE; ++j){
			dTmp += gMatrixA[i][j]*gVectorX[j];
		}
		dTmp-=gVectorB[i];
		dL2Norm += pow(dTmp, 2.0);
	}	
	*pL2Norm = sqrt(dL2Norm);

	return 0;
	
}
/********************************************/
void showMatrix()
{
	int i, j;
	for(i=0; i<SIZE; ++i){
		for(j=0; j<SIZE; ++j)
			printf("%f\t", gMatrixA[i][j]);
		printf("\n");
	}
	printf("\n");
	for(i=0;i<SIZE; ++i){
		printf("%f\n", gVectorB[i]);
	}
	printf("\n");
	
	for(i=0;i<SIZE; ++i){
		printf("%f\n", gVectorX[i]);
	}
	printf("\n");


}

int main()
{
	double dL2Norm=0.0;

	init();
#if 0
	showMatrix();

	gaussian_elim();
	backSubstitution();
	getL2Norm(&dL2Norm);
	printf("L2Norm: %.12f\n", dL2Norm);

	showMatrix();
#endif 
	fini();
	return 0;
}
